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Airborne wind energy systems are capable of extracting energy from higher wind speeds at higher altitudes. The configu¬ 
ration considered in this paper is based on a tethered kite flown in a pumping orbit. This pumping cycle generates energy 
by winching out at high tether forces and driving a generator while flying figures-of-eight, or lemniscates, as crosswind 
pattern. Then, the tether is reeled in while keeping the kite at a neutral position, thus leaving a net amount of generated 
energy. In order to achieve an economic operation, optimization of pumping cycles is of great interest. 

In this paper, first the principles of airhome wind energy will he briefly revisited. The first contribution is a singularity- 
free model for the tethered kite dynamics in quaternion representation, where the model is derived from first principles. 
The second contribution is an optimal control formulation and numerical results for complete pumping cycles. Based 
on the developed model, the setup of the optimal control problem (OCP) is described in detail along with its numerical 
solution based on the direct multiple shooting method in the CasADi optimization environment. Optimization results for 
a pumping cycle consisting of six lemniscates show that the approach is capable to find an optimal orbit in a few minutes 
of computation time. Eor this optimal orbit, the power output is increased by a factor of two compared to a sophisticated 
initial guess for the considered test scenario. 
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1 Introduction 

The idea of airborne wind energy (AWE) is to generate usable power by airborne devices. In contrast to towered wind 
turbines, airborne wind energy systems are flying, usually connected by a tether to the ground, like kites or tethered 
balloons. AWE systems exploit the relative velocity between the airmass and the ground and maintain a tension in the 
tether. They can be connected to a stationary ground station, or to another moving, but non-flying object, like a land or sea 
vehicle. Power is generated in form of a traction force, e.g. to a moving vehicle, or in form of electricity. The three major 
reasons why people are interested in airborne wind energy for electricity production are the following; 

• Eirst, like solar, wind power is one of the few renewable energy resources that is in principle large enough to satisfy 
all of humanity’s energy needs. 

• Second, in contrast to ground-based wind turbines, airborne wind energy devices might be able to reach higher al¬ 
titudes, tapping into a large and so far unused wind power resource Q. The winds in higher altitudes are typically 
stronger and more consistent than those close to the ground, both on- and off-shore. 

• Third, and most important, airborne wind energy systems might need less material investment per unit of usable power 
than most other renewable energy sources. This high power-to-mass ratio promises to make large scale deployment of 
the technology possible at comparably low costs. 

In order to achieve an economic and competitive operation of airborne wind energy (AWE) setups, optimal control of 
energy production is the key to success. In this paper, a novel formulation of the dynamics of a specific AWE system is 
given, and it is shown how it allows one to solve challenging optimal control problems for this system. 

• Corresponding author, e-mail: michael.erhard@skysails.de, Phone +494070299 203, Fax +494070299 222 
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The paper is organized as follows: in Section]^ we introduce the main ideas of airborne wind energy, following the lines 
of 0, and then explain in Sectionj^the experimental system developed by the company SkySails, the modeling of which 
is the main subject of this paper. Starting point is an existing differential equation model of the system from Q, described 
in Section which suffers from singularities caused by the choice of coordinate system. The main contribution of this 
paper is presented in Section]^ where a new singularity free model based on quaternions is developed. It is demonstrated 
by numerical simulation that the old and new model describe identical system dynamics when away from the singularities, 
while only the new model can be reliably simulated everywhere. It is then shown in Section|^how the new quaternion-based 
model can be used to formulate the complex periodic optimal control problem (OCP). In Sect.|^ the discrete formulation 
of the OCP to be solved in the dynamic optimization environment CasADi Q is given. Furthermore, resulting complete 
operational pumping cycles for the SkySails airborne wind energy system are presented. The paper ends with a conclusion 
in Section in 


2 Crosswind kite power and pumping cycles 

Every hobby kite pilot or kite surfer knows this observation: As soon as a kite is flying fast loops in a crosswind direction 
the tension in the lines increases significantly. The hobby kite pilots have to compensate the tension strongly with their 
hands while kite surfers make use of the enormous crosswind power to achieve high speeds and perform spectacular stunts. 
The reason for this observation is that the aerodynamic lift force Fj, of an airfoil increases with the square of the flight 
velocity, or more exactly, with the apparent airspeed at the wing, which we denote by Va,. More specifically, 

= ]^pACi,vl, (1) 

where p is the density of the air, A the airfoil area, and Cl the lift coefficient which depends on the geometry of the airfoil. 

Thus, if we fly a kite in crosswind direction with a velocity Va that is ten times faster than the wind speed the tension 
in the line will increase by a factor of hundred in comparison to a kite that is kept at a static position in the sky. The key 
observation is now that the high speed of the kite can be maintained by the ambient wind flow, and that either the high 
speed itself or the tether tension can be made useful for harvesting a part of the enormous amount of power that the moving 
wing can potentially extract from the wind held. 

The idea of power generation with tethered airfoils flying fast in a crosswind direction was already in the 1970s and 
1980s investigated in detail by the American engineer Miles Loyd | fT8| . He was arguably the first to compute the power 
generation potential of fast flying tethered wings - a principle that he termed crosswind kite power. Loyd investigated (and 
also patented) the following idea: an airplane, or kite, is flying on circular trajectories in the sky while being connected 
to the ground with a strong tether. He described two different ways to make this highly concentrated form of wind power 
useful for human needs, that he termed lift mode and drag mode: while the lift mode uses the tension in the line to pull a 
load on the ground, the drag mode uses the high apparent airspeed to drive a small wind turbine on the wing. 

This paper is concerned with systems operating in lift mode, which has the advantage over the drag mode that that it 
does not need high voltage electrical power transmission via the tether. The idea is to directly use the strong tether tension 
to um'oll the tether from a drum, and the rotating drum drives an electric generator. As both the drum and generator can 
be placed on the ground, we call this concept also ground-based generation or traction power generation. Lor continuous 
operation, one has to periodically retract the tether. One does so by changing the flight pattern to one that produces much 
less lifting force. This allows one to reel in the tether with a much lower energy investment than what was gained in the 
power production phase. The power production phase is also called reel-out phase, and the retraction phase reel-in or return 
phase. Loyd coined the term lift mode because one uses mainly the lifting force of the wing. But due to the periodic reel-in 
and reel-out motion of the tether, this way of ground-based power generation is often also called pumping mode', sometimes 
even the term Yo-Yo mode was used to describe it. 

It is interesting to compare crosswind kite power systems with a conventional wind turbine, as done in Lig. [T] which 
shows a conventional wind turbine superposed with an airborne wind energy system. Seen from this perspective, the idea 
of AWE is to only build the fastest moving part of a large wind turbine, the tips of the rotor blades, and to replace them by 
automatically controlled fast flying tethered kites. The motivation for this is the fact that the outer 30% of the blades of a 
conventional wind turbine provide more than half of the total power, while they are much thinner and lighter than the inner 
parts of the blades. Roughly speaking, the idea of airborne wind energy systems is to replace the inner parts of the blades, 
as well as the tower, by a tether. 

The power P that can be generated with a tethered airfoil operated either in drag or in lift mode had under idealized 
assumptions been estimated by Loyd | fT8| to be approximately given by 
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Fig. 1 The basic idea of the power 
generating device is to only build the 
wing tips of a gigantic windmill, in 
form of tethered airfoils, i.e. kites. Fast 
crosswind motion of the kites creates an 
enormous traction force [graphics cour¬ 
tesy Boris Houska], 



Fig. 2 SkySails functional prototype setup for kites of sizes ranging 
from 20 to 40 m^ (30 m^ shown here). A tether in the range 150-400 m 
transfers the forces from the flying system to the ground based main 
winch, which is attached to a motor/generator. A specific feature of this 
setup is the single tether and the airborne control pod located under the 
kite. The steering actuator in the control pod pulls certain lines in order 
to steer the kite. 


where A is the area of the wing, Cl the lift and Cd the drag coefficients, and the wind speed. Note that the lift-to-drag 
ratio E = ^ enters the formula quadratically and is thus an important wing property for crosswind AWE systems. For 
airplanes, this ratio is also referred to as the gliding number, it describes how much faster a glider without propulsion can 
move horizontally compared to its vertical sink rate. 

Since the first proposal of using tethered wings for energy harvesting by Loyd, academic research and industrial de¬ 
velopment activities in the field of airborne wind energy (AWE) have been significantly increasing, especially during the 
last 10 years. The main advantage of the AWE technology is that airborne systems are capable of harvesting energy from 
higher wind speeds at higher altitudes. In addition, as most setups require less installation efforts compared to conventional 
wind turbines, e.g. foundations and towers, this technology is a promising candidate to become an additional source of 
renewable energy in the near future. For an overview on the different geometries which have been developed so far, the 
reader is referred to (Tg and the monograph on AWE |[T). 

3 The SkySails Power prototype 

The geometry under consideration in this paper is based on the SkySails Power prototype, shown in Fig. This setup 
features an airborne control pod located under the soft kite connected to a ground based winch with motor/generator by 
a single tether. The functional prototype successfully demonstrated fully automated energy generation by autonomously 
controlled pumping cycles GD- The principle of energy generation is depicted in Fig. and shall be briefly summarized; 
during the power phase, the kite is flown dynamically in lemniscates (figures of eight). This dynamical crosswind flight 
leads to high tether forces. By reeling out the tether, electrical energy is produced. At a certain line length, the kite is 
transferred to the neutral zenith position. At this low-force position, the tether is reeled in (return phase) consuming a 
certain portion of the previously generated energy. Finally, a reasonable amount of generated net power remains. Hence, 
the average power output over complete cycles will be the target of the optimization implementation presented in this paper. 

The optimal control is based on a dynamical model of the system, and the chosen numerical techniques are similar to 
the ones used in fTS) . In order to allow for optimization of complete pumping cycles, a simple model should be chosen in 
order to keep the number of optimization variables low. The model developed in this paper is based on a very simple model 
set up for the SkySails system Q, which describes the basic dynamic properties of the system and has been experimentally 
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2. Transfer Phase 


1. Power 
Phase 



3. Return 
Phase 


Fig. 3 Trajectory of an experimentally flown pumping cycle (TT) 



Fig. 4 Coordinate system definitions of the fixed reference vectors ex,&y,&z and the body 
frame vectors eroii, Spitch and Oyaw The position is determined by the two angles 93 , i9 and 
the tether length 1. The wind direction is defined in e^^-direction. 


validated |10|. The quaternion representation allows for singularity-free equations of motion, and the replacement of 
trigonometric functions seems to be advantageous for the optimization process. 

A further advantage of the considered system is that an experimentally flown trajectory (see e.g. Fig. can be used as 
initial guess and the optimization result can be directly compared to it in return. It should be noted that the optimization has 
to be subject to certain geometrical constraints, which will be discussed in this paper. Particularly, topological constraints 
have to be imposed to preserve the topology of the pumping cycle during the optimization process (the 6 lemniscates given 
in the initial guess, compare Fig.[^. 


4 Setup and reference kite model 

In this section, the model from ®,|Tg is briefly summarized. The coordinate system is defined in Fig. The state vector 
of the dynamical system consists of the 3-dimensional orientation of the flying system, given by {tjj, ip, z?} and the tether 
length 1. 

x= (3) 


Note that the tether of the system leads to a relation between orientation and position, which is given in the basis { e^;, Oj,, } 

by: 


r = I 


cos'd 

sin ip sin id 
— COS ip sin id 


(4) 


The wind direction is defined in e 3 ;-direction. In order to provide a clear picture of the coordinate system definition, a 
comparison with the navigation task on earth shall be discussed as illustrated in Fig.|^ The coordinate system corresponds 
to an earth coordinate system with symmetry (rotational) axis in wind direction, which could be achieved by rotating the 
north pole by 90 degrees towards the wind direction. According to the angle mapping given in Table [T] the two angles 
ip, -d correspond to longitude and latitude, hence determine the position on earth. The orientation angle ij} corresponds to 
bearing, i.e. direction w.r.t. north in our earth system. This angle determines the direction of motion (course) and the angle 
of ■!/) = 0 describes heading directly against the wind. It should be Anally noted that due to the ambient wind, the angle 
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Fig. 5 Demonstrative explanation of the used Euler angles 
by comparing them to angles used for navigation on earth. 
For that purpose, the earth’s rotation axis has to be tilted by 
90 degrees with the north pole towards the wind direction. 
The angles ip, correspond to the position on earth and t/) 
to the bearing. The exact mapping is given in Table [T] 


Table 1 Mapping of the Euler angles of the tethered system to 
earth navigation quantities (compare Fig.[^. 


angle in model 

’earth’ navigation quantity 


earth longitude 

t? 

earth latitude 

V' 

bearing (direction w.r.t. north pole) 


does not exactly coincide the kinematic flight direction (course). This effect is fully covered by the developed model and 
further explained in mi- 

The input vector of the model comprises the inputs of the two control actuators: the steering deflection of the airborne 
control pod, 5 and the winch speed of the tether respectively. 


u = 




^(winch) 


T 


(5) 


In the following, the equations of motion shall be summarized for reference purpose only. This model was derived 
based on the same assumptions as those that will be used in the subsequent section for the derivation of the quaternion 
representation. As conducted in detail in Q, pO) , the equations of motion read: 


-0 

= gkVaS + ip COS l} 

(6) 


Va . 

= ; . .sinV^ 

/ sinv 

(7) 

d 

. Q 1 1 

= —— Sin V 

(8) 

i 

— .j; (winch) 

(9) 

Here, the airpath speed of the kite, Va, is given by 


Va 

= v^E cosd — IE. 

(10) 


The set of equations involves three parameters: lift-to-drag (glide) ratio E, steering response proportionality constant 
and ambient wind speed Uw- Further explanations of parameters including values used later in the paper are summarized in 
Table 1^ Finally, the tether force can be computed based on the air path speed Ua by 


Ftether 


qACs. E 2 
—^- / 

2 ^/lEE^^ 


( 11 ) 


with density of air q, projected kite area A and aerodynamic force coefficient Cr. 


5 Quaternion based kite model 

In the following, the equations of motion shall be derived in a singularity-free formulation based on quaternions. Quater¬ 
nions are a well-known tool in aerospace modeling and estimation to avoid the singularities of Euler angles fTT] . An 
alternative singularity-free representation could be based on rotation matrices as natural coordinates directly, as is shown 
for modeling of AWE systems in El- Note that the rotation matrix in the following will be used only as intermediate step 
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Table 2 Description of system parameters and values used for optimization 


symbol 

value 


description 

A 

21.0 

m^ 

projected area of the kite 

Cr 

1.0 


aerodynamic force coefficient, given by Cr, = 
\JC'rE Cp for aerodynamic lift (drag) coefficient Cl 
(Cd) 

^max 

0.6 

1/s 

maximal control actuator steering speed 

^max 

0.7 


maximal control actuator steering deflection 

E 

5.0 


lift-to-drag (glide) ratio. Note, that in contrast to most 
other systems, the SkySails system is operated at constant 
glide ratio 

5k 

0.1 

rad/m 

proportionality constant relating turn rate of the kite to 
steering actuator deflection 

^max 

300 

m 

maximal available tether length 

^niin 

0.35 

rad 

minimal elevation angle w.r.t. the horizontal plane 

Q 

1.2 

kg/m^ 

air density 

^a,min 

5.0 

m/s 

minimal air path speed in order to guarantee flight sta¬ 
bility and keep system tethered (this corresponds to free 
flight velocity) 

(winch) 

V ■ 
min 

-5.0 

m/s 

typically chosen as —Uw/S, leads to a windward limit of 
« 110° 

nw 

10 

m/s 

ambient wind speed 


Reference frame axes 


eroll 


Wind 


^pitch 


[Kite body frame axes ] 


Fig. 6 Reference frame axes for the derivation of the quaternion 
equations of motion. The depicted orientation corresponds to the 
zero-rotation, i.e. i? = / or q = [1,0,0,0]. 


in the derivation towards quaternions. The base coordinate system is depicted in Fig. Si. The time dependence w.r.t. to this 
reference frame is expressed by rotation matrices R{t) G SO(3), i.e. 


eroll(i) = 

-R{t)Gz 

(12) 

6pitch(^) — 

—R{t)ey 

(13) 

0 

II 

-R{t)ex 

(14) 


The zero-rotation state is defined by i? = /, corresponding to Oyaw = —Gx, Gpitch = —Gy and OroU = —Gz as illustrated 
in Fig.|^ Note that this seems to be an arbitrary choice a priori. However defining the ambient wind direction, given by e^, 
in the kite’s yaw-axis leads to a beneficial simplicity of the equations of motion later. 

It should be emphasized that due to the nature of the tethering, position and orientation are related. Hence, the position 
of the kite is given by 


r(t) = -i(f)eyaw(f) = l(t)R(t)ex 


(15) 


The air path speed vector can be calculated as sum of the ambient wind vector Vv,Gx and reversed kinematic velocity vector 
—r (apparent wind opposed to motion). Using the time derivative of (15 1 and introducing the speeds Uroii and Wpitch in the 
tangent plane, compare Fig. yields 


i(f) — '^wGx npitch(f)^pitch(f) ^roll(f)tlroll(f) T Kf)®yaw(f) 


(16) 
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top view 


^pitch 


®roll 





V 


a 



side view 



Fig. 7 The kite is assumed to always fly in an aerodynamic equilibrium 
state. As aerodynamic model, absence of transversal air flow (top view) and 
constant lift-to-drag-ratio (side view) is demanded. Figure taken from ||10| 


Fig. 8 Motion in the tangent plane. Due to the tether, the velocities in the tangent plane Uroii and Upitch 
determine Wroii and ujpitch as given in \22) and i[2'i). 


5.1 Model assumptions 

For AWE systems, the aerodynamic forces are typically large compared to masses, which in fact is a prerequisite for 
tethered flight. This is particularly true for the SkySails Power system. Therefore, acceleration effects play a minor role 
and can be neglected. In addition, the system is assumed to fly always in its aerodynamic equilibrium state. These two 
simplifications lead to the simple structure of the equations of motion. 


5.1.1 Aerodynamics 

The aerodynamics of the system is reduced to two conditions, depicted in Fig. and given in the following. 

1. The kite experiences no transversal wind flow, i.e. no ’side-slip’ angle occurs, compare Fig.|^. 

(epitch(f), Va(f)) = 0 (17) 

Inserting ( [Thl l yields 

(^pitch(f )5 '^pitch(f) ~F ^ (^pitch (^)5 ®yaw(^)) — 0 (18) 

-V-^ 

=0 

Thus, the condition results in 

^pitch(f) (19) 


2. The lift-to-drag ratio is constant and given by the parameter E. Hence, the kite is assumed to fly always at the same 
angle of attack. Considering the geometry of Fig. |^) yields 


(erou(t), Va(f)) = E{e yaw W>Va(i)) 

Inserting ( [T6] l produces 

^roll(^) — '^w (7?(f)(F'Oa; Eli^t) 


( 20 ) 


( 21 ) 


Finally, by considering the motion in the tangent plane, compare Fig. the following relations for the turn rates in kite- 
body-frame can be derived. 

^pitch(f) 


Wroll = 


= “7777 


l{t) l{t) 


^roll(0 /D/j.\/77i \ \ I 771 

t^pitch - --j^{Rit){Ee,-e,),e,)+E — 


( 22 ) 

(23) 
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5.1.2 Steering 

Steering of the kite is accomplished by the control actuator in the control pod. Applying an actuator deflection 5 to the 
steering lines tilts the kite canopy resulting in a curve flight or turn rate around the Oyaw-axis. The steering behavior can be 
described by the following equation 

Wyaw = g\-iVj ( 24 ) 

This turn-rate law states that the rotation rate is proportional to a constant system parameter gy, the deflection 5 and the 
scalar air path speed Ua, which is defined as component of Va in —ej-oirdirection 

= - (eroii(i), Va(f)) = {R{t)e^, ( 25 ) 

The validity of the turn-rate law has been experimentally shown for different kites In addition, a nice 

derivation based on geometric considerations can be found in GD- It should be mentioned that correction terms for mass 
effects 0 and aerodynamic effects might be added, but are neglected in this paper. 


5.2 Quaternions 

The final formulation of the equations of motion is based on quaternions fTT), introduced as 


q = 


90 

91 

92 

93 


( 26 ) 


The time evolution as function of turn rates, which are given in the body frame, reads 



0 

tva; 

-UJy 

-UJz 

1 

tv a, 

0 

UJz 

-UJy 

2 

Uy 

-UJz 

0 

tv a: 


_ 

UJy 


0 


(27) 


Mapping of the kite axes to uix,ujy, oj^ has to be considered for R = I. Regarding (12 14 1 or Fig. yields ojx = — 


UJy 


UJ, 


, = —Wpitch and ojz = — Wroii- The quaternion time evolution then results in 



0 

^yaw 

^pitch 

Unroll 

^yaw 

0 

Unroll 

^pitch 

^pitch 

^roll 

0 

^yaw 

^roll 

^pitch 

^yaw 

0 


The relation to the rotation can be given by 



■ 9o + 9? - 9i - 9| 

2(9192 - 9093 ) 

2(9193 + 9092 ) 

R{t) = 

2(9192 + 9093 ) 

9o - 9? + 9i - 9 I 

2(9293 - 9o9i) 


. 2(9193 - 9092 ) 

2(9293 + 9o9i) 

9o - 9? - 9 I + 9i - 


(28) 


(29) 


The derivation of the system’s equations of motion is done as follows: The turn rates in kite-body-frame, Wroih Wpitch and 
Wyaw are expressed as functions of q and subsequently inserted in ( [2^ . In order to achieve this, (29 1 is used for R(t) and 
inserted into ( |22] i, ( |2^ and ( |24| i. 


5.3 Equations of motion 

The state vector for the quaternion formulation is given by 

[ 90 , 91 , 92 , 93 ,^]^ 

and the input vector by 

5,u(winch)l^ 


( 30 ) 


(31) 
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Performing the derivation steps of the previous section yields the following set of equations of motion: 


-92 


90 ( 9 ! + 9i) 


9i 

-93 

9o 

1 

91 ( 9 ! + 9I) 
- 92(90 + 9?) 

2 

1 1 

OJ 0 

9i 


- 93(90 + 9 ?) 


92 


i = 

with air path speed 

Ua = Ev^{ql + ql-ql- qj) - El 


(32) 

(33) 

(34) 


It should be remarked that if the winch speed is defined relative to the wind speed by substituting I = the 

product of wind speed Uw and time t (or cycle time T) becomes an invariant. As a consequence, solutions for different 
wind speeds are simply obtained by scaling the dynamic’s time axis proportional to 1/v^. 

Although these equations of motion (EQM) preserve the quaternion norm (d/dt)||q(f)|| = 0, a stabilization term in 
order to compensate for numerical inaccuracies should be added. The EQM for the numerical optimization read then 


q = <!>(■)-7,(l|q|P-l)q 


(35) 


where $(•) denotes the equations of motion ([^ 34i and the damping constant jg is chosen to lead to a slow decay (jg = 
0.011/s for the parameters given in Table|^. 


5.4 Transformation to Euler angles for comparison 

Eor sake of completeness and in order to allow for a simple geometric interpretation of the results in the Cartesian coordinate 
system (see Eig. |^, transformations from the set of Euler angles {(p, i}, tp} to quaternions and back shall be summarized. 
The sequence of rotations 


R = R,{p)Ry{d)R,{-P;) 

is considered in quaternion representation 



■ cos (f) ■ 


■ cos(f) ■ 


cos ) 

q = 

sin(f) 

0 

0 

0 

sin(f) 

0 

-sin(l) 

0 


0 


0 







0 


Calculating the quaternion product (0) 113 of the single rotations yields 

cos (f) cos (I) cos (l) + sin (f) cos (f) sin 
sin (f) cos (f) cos (f) - cos (f) cos (f) sin (f J 
- sin (f) sin (f) sin + cos (f) sin (f) cos 
sin (f) sin (f) cos + cos (f) sin (f) sin 

The kite position w.r.t. {e^:, 0 ^, 6 ^} can be found by inserting into as 


r = I 


(9o + 9? - - 9i) 

2(9093 + 9192) 
2(9193 - 9092) 


The Euler angles can be computed from the quaternions by 

p = atan2 ((^093 + 9192 ), (9092 - 9193 )) 
■d = arccos(g^ + 9? - 9^ - 9^) 

Ip = atan2 ((go93 - 9192 ), (9092 + 9193 )) 


(36) 


(37) 


(38) 


(39) 


(40) 

(41) 

(42) 


where atan2 is the four-quadrant arctangent function. 
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Fig. 9 Simulation of the simple model ODE, formulated in Euler angles (solid red) and quaternions (dashed blue). The 
3d trajectory is shown on the left, the right plot zooms into the vicinity of the singularity, where the Euler angle formulation 
fails. The points in the right subfigure represent the discrete simulation timesteps of constant duration. 


5.5 Simulation example 

In order to complete the presentation of the quaternion-based model, a simulation example will be discussed. The simula¬ 
tions have been carried out using the quaternion-based model on the one and the reference model from Sect. |^on the other 
hand. Both computations have employed a Runge-Kutta-4 ODE propagation using timesteps of r = 0.1 s, a tether length of 
I = 100 m and the initial conditions (p(0) = 0, t?(0) = arctan(i?) and ipi^) = 0. A constant steering command of <5 = 0.186 
has been used to fly next to the singularity of the Euler angle representation. The simulation results are depicted in Eig. 

It can be recognized that the propagation of the model (|6]jT0|) fails due to the singularity, while the quaternion formulation 
([32||34|) leads to correct numerical simulation results. 


6 Optimal control problem formulation 

The optimization task is to control power cycles with optimal power output averaged over complete cycles. This quantity 
is computed as integral over the mechanical power 

1 /■; , E If;,, 

^ ^tether df- df (43) 

0 0 

It should be mentioned at this point that an open loop optimal control problem is set up, which is based on controls as 
model input, but does not involve any feedback control based on measured values. 


6.1 Setup of the optimization problem 

As steering speed will also be limited, the steering deflection 5 is introduced as additional state to the model and the pod 
steering actuator speed, defined as 5s, is used as control, instead. Eurther, the state W with W{Q) = 0 for computing the 
integral in (43 1 is added. Thus, the following relations have to be added to complete the set of EQM ([^ 34 1 


5 = 
W = 


5s 

iv^ 


(44) 

(45) 


The control vector now reads 


u= 5s, V 


(winch) 


(46) 


and the state 


X = [W,5,l,qo,qi,g2,q3]'^ 


(47) 
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The NLP problem is formulated as minimization of 


/ = -^ 


iv^ df + ea / dt = — 


W{T) 


es / |(5(f)p dt 


(48) 


0 0 0 

The second term introduces a penalty in order to achieve a smooth steering behavior. The weighting factor es has to be 
chosen appropriately. 

The optimization, i.e. minimization of /, has to be carried out by variation of x(f),u(f),r subject to the following 
constraints. First the optimization is subject to the model equations of motion (32 -[34|. In addition, physical, geometrical 
and topological constraints have to be added, which will be discussed in the subsequent sections. 


6.2 Physical constraints 


Due to weight and geometrical design considerations for the airborne system, the two limitations of limited steering speed 
and deflection of the control pod actuator have to be taken into account. Hence, the steering speed limit is given by 

l^sl < (49) 

and the limited steering deflection is implemented by 

1^1 < <5max (50) 


respectively. In addition, winching speed is limited 

(winch) (winch) 

‘^min — ^ 


(51) 


Note that this bound on the reel-in speed imposes implicitly an constraint on the windward position during the reel-in phase, 
as can be shown by computing the steady state angle lTg,i[n|. Further, in order to keep the system tethered, a certain 
minimum air path speed and hence tether force is required. Recalling ( [T0| yields 

tia — d~ Qi f ?2 ^ 3 ) — tia.min (52) 


Note that this condition is needed to guarantee the validity of the model during the optimization process. 


6.3 Geometric constraints 

As already introduced, optimization is done for closed pumping orbits, which impose periodic boundary conditions to the 
optimization problem 


<l{T) 

= q(0) 

(53) 

1{T) 

= m 

(54) 

6iT) 

= S{0) 

(55) 


We also add one non-periodic boundary condition in order to define the initial value of the “energy state” W: 

VF(0) = 0 (56) 

The tether length has to be constrained in order to keep the cycle within a certain range of line length. This is accomplished 
by setting an upper bound to the state 

I < ?max (57) 

For a real world system, earth (or water) surfaces are a serious constraint on the equations of motion. Practical safety 
considerations taking into account flight trajectory deviations due to wind gusts recommend the choice of a certain minimum 
elevation angle 0niin- Geometric considerations yield 

>tan0„,i„ (58) 

Using ( |42| i results in the constraint 

(9o + 7i - 92 - 93 ) tan0i„in + 2 ( 91(73 - 9092 ) < 0 (59) 

It should be finally noted that the operational flight altitude strongly depends on the wind profile and ideally, the optimal 
flight altitude lies above the given safety limit. However, for the optimal control computations in this paper, the wind held 
is assumed constant and homogeneous for simplicity and therefore these constraints play a major role. 
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Fig. 10 Topological constraints shown for a pumping cycle consisting of two lemniscates. In order to avoid ’unwinding 
of the topology’ by the optimizer, the trajectory horizon is divided into stages and additional constraints are imposed on 
those to preserve the topology. 


6.4 Topological constraints 


A crucial point in the formulation of the optimization problem is the topology of the trajectory. For operational reasons, 
basically in order to avoid twisting of the tether, the patterns of choice are lemniscates. Hence, one wishes to preserve the 
topology of the pattern during the optimization process by imposing certain constraints. If this is not done, the optimizer 
tries to ’unwind’ certain trajectory features, which in almost all cases leads to a failure of the optimization process. 

The topological constraints aim at keeping the configuration of the pattern, e.g. two lemniscates as sketched in Fig.fTO] 
or six, as in the actual computations of this paper. The trajectory can be divided into stages, mainly parts featuring the 
property flying to the right (left), which corresponds to (^> 0 and ip<Q, respectively. Formally, an even number of N time 
intervals with duration Ti is introduced. Further, the time points to,... ,t]s[ with the property Ti = ti — ti-i are defined as 


shown in Fig. 10 Note that fo = 0 and that = T is the overall pumping cycle time. 

As (p is not a state of the system, 0 can be used to formulate an equivalent constraint on ijj. Considering the sign of 
sin Ip together with the behavior of the four-quadrant arctangent function atan2 in (40 1 , the following conditions can be 
deduced 


<?O93-9i 92>0 t(2i) < ^ < i(2i+l) * = 

9093 - 9192 < 0 f(2i_l) < t < t(2i) 

The upper (lower) relation corresponds to {<p < 0). 


6.5 Summary of the optimal control problem 

The optimal control problem (OCP) can be stated formally as follows: 


minimize 


tN 

E{x{tN),tN) + J L{u{t)) dt 

0 


( 61 ) 


subject to 

x(<)-$(x(f),u(f)) = 0, 

rboundary(^(fo)5 ~ 

< 0 , 

hi{x{t)) < 0, 


t G [fo,fAr], 
t G 

t G i = l,...,N. 


(ODE model) 

(boundary conditions) 

(path constraints) 

(stage wise path constraints) 
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• integration steps 



K control steps j t in 


I multiple shooting intervals 


rii controls u; j and states Si j 


• stages 


+ 


+ 



Fig. 11 Setup of the discretization grid. 


N interval times Ti 

—I- 


with 


E{x{tN),tN) 


WjtN) 

In 


and L{\i{t)) = 


(62) 


The set of equations for the ODE model $ (•) is defined by (28 34 1 and (44 451. The boundary conditions ^boundary (•) ^6 
given by the periodicity and initial state constraints (53 -56 1 , the path constraints h{-) by (49 52i, (57 1 , (591 and the stage 
wise path constraints /li(-) by (60l. 


7 Numerical computation of optimal pumping orbits 

In order to numerically solve the optimal control problem, a direct method is used to transfer the continuous time OCP into 
a nonlinear programming problem (NLP). This is accomplished by using the direct multiple shooting method. The direct 
multiple shooting method that was originally developed by Bock and Plitt Q performs first a division of the time horizon 
into subintervals, that we will call multiple shooting (MS) intervals in the sequel, and it uses piecewise control discretization 
on these intervals. The continuous time dynamic system is transformed to a discrete time system by using an embedded 
ordinary differential equation (ODE) solver to solve the ODE on each MS interval individually. Details of the setup of 
the discrete time control problem will be summarized in the following subsections, and finally numerical results will be 
presented, that have been obtained by using the optimization environment CasADi @13 and the NLP solver IPOPT 


7.1 Multiple shooting and control discretization grid 

The choice of time discretization grid size choice is a compromise between keeping the total number of optimization 
variables low and at the same time reproducing the dynamics of the controls and states in an adequate way. In order 
to represent controls and states of different timescales, nested grids are introduced as depicted in Pig. 11 The upper 


level grid is the splitting of the trajectory into stages as introduced in Sect. |6.4| The stage of duration Ti is divided into Ui 
equidistant multiple shooting (MS) intervals for the states and the winch control. Note that the stages can comprise different 
numbers rii of MS intervals in order to account for different phases as power and return phase, respectively. These MS 
intervals are divided into K subintervals for the kite steering actuator speed (control) to allow for a proper consideration 
of the constraints ( |49] l in the solution. The initial value on teach the control vector on each MS interval consists of iT + 1 
components, i.e. 




(winch) 

■ 

^,3 


(63) 


The idea of multiple shooting is to solve the ODE on small intervals with duration starting with initial values Si j The 
used scheme is sketched in Pig. 12 The integration is performed from an initial state s, j to the final state x* j, using 
the controls j . Here, we are utilizing K integration steps of the classical Runge-Kutta integrator of order 4 (RK4), each 
with a step size ti = In order to give a formal representation, the single RK4 integration step $hk : s —>^ x shall be 
defined as 


(64) 


x = 4>rk(a,u(“),<5,s) 
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Fig. 12 Principle of the implemented shooting scheme for K — 3. The ODE model is solved on intervals by K integration 
steps, starting at Sij and ending at Xij. The trajectory is ’connected’ by imposing constraints, e.g. S; (j+ij — Xij = 0. 
Note that in contrast to the usual implementation of multiple shooting, several piecewise constant control steps are taken 
per MS interval for one of the controls. 


The result of the K RK4 integration steps on the interval with duration — is then given by 


^i,jiSi,j,Uij) = $rk(Tj, , ^'RK(rj, 77 ”™^ , ^'RK(rj, 7 . 7 ”™^ Si,i) • • • )) 


(winch) h 


(winch) h 


K integration steps 


(65) 


It should be emphasized that separate piecewise constant control steps on the subgrid are provided for the integration 
steps of the MS interval. 

In summary, the optimization vector w = [w7..., can be subdivided into subvectors with all variables 

belonging to stage i, which are given by 




iT 


( 66 ) 


7.2 Summary of the discretized optimal control problem 

The complete discretized OCP can be stated as follows 


minimize £^(sAr,n„, Ti) + T,) 


(67) 


subject to 

X 



i = 1,. 

..,N;j = 0,. 

. .,(n*-l) 

(continuity between MS intervals) 

Si.ri; — S(i_|_i)_0 = 0 

i = 1,. 



(continuity between stages) 

^boundary (®1,0) ®A^,niv ) “ 6 




(boundary conditions) 

h{sij,Uij) < 0 

i = 1,. 

..,N;j = 0,. 

. .,(n*-l) 

(path constraints) 

< 0 

i = 1,. 

..,N;j = 0,. 

. .,(n*-l) 

(stagewise path constraints) 


with 




K 

sr 




T, 




winch) (winch) 
3 '^next(2,j'' 


( 68 ) 
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where next(i, j) denotes the subsequent index pair in time. Note that the second term penalizes changes in winch speed 
(i.e. acceleration) and has been added to the discretized OCP only. The discretized path constraints h{-) evaluate the 
constraints (49 52 1 , (57 1 , (59 1 on the multiple shooting grid, and the discretized stage wise path constraints hi{-) evaluate 
(|60)l on the grid. 


7.3 Numerical results 

In the following, numerical results for optimal complete pumping cycle will be presented. The parameters are given in Table 
l^and are related to the SkySails functional prototype of Fig. As initial guess, an experimentally flown trajectory is used. 
The topology of 6 lemniscates leads to = 12 and the (initial) cycle time of T « 170 s is divided up into — 250 

intervals. The number of RK4 integration steps is chosen as AT = 3 on each MS integration interval. Thus, the optimization 
vector contains a total of 2762 optimization variables. The 3d trajectories at different iteration steps are shown in Fig. 13 
It can observed that optimization starts far from optimum and runs through some crude path configurations towards the 
optimal solution. This behavior could be attributed to the interior point method of the used IPOPT solver. 

The optimization figures of merit are given by the Loyd factor pLoyd, which is defined by 


P 


^Loyd — 


Pt 


Loyd 


(69) 


This factor corresponds to the Loyd limit, which corresponds to the generated power Puoyd for a kite flying ’maximal 
crosswind’ (at rl = 0) continuously, i.e. no retraction phase, compare © and fTS) . For the derived model, this quantity is 
given by 


PLoyd — 


pCrA 4p2 E 
2 27 v/T+A^' 


(70) 


The according time series for the controls and power, given as Ptether^, are presented in Fig. 


14 


In summary, it can be stated that the optimizer increased the power generation of the pumping cycle by more than a factor 
of two from pLoyd « 0.15 to PLoyd « 0.33. This is partly due to the fact that the optimized figure lies deeper in the crosswind 
region. Further, transfer and return phases are carried out more efficiently. The lemniscates are deformed to minimal curves 
fulfilling the imposed constraints. For sake of a fair assessment, it should be noted that the experimental trajectory was 
taken from a test flight aiming at flight controller development rather than optimal power output. Nevertheless, valuable 
suggestions for improvement of even these ’classical’ control setups can be extracted from such optimization results, see 
e.g. the winch controller implementation for the transfer phase in mi- 


8 Conclusion 

In conclusion, a singularity free model for tethered kite dynamics based on quaternions has been derived, which can 
be broadly applied for simulation and optimization purposes. It has been demonstrated that optimization runs of complete 
pumping cycles within a few minutes are feasible with this model. This opens up the way to further extended and systematic 
studies of AWE efficiencies. It should be finally remarked that these might most likely involve further extensions to the 
model to take into account e.g. mass effects, tether drag and a wind profile. 
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